randn('state',32);
rand('state',32);
load c:/data/restat/analysis_data/postoillogitdata;
id=data(:,1);
s_jt=data(:,2);

x1=[data(:,3:36)];
clear data;
load c:/data/restat/analysis_data/postoilblpiv;
iv=data(:,3:28);
brand=floor(id/100000);
company=floor(brand/10000);

pre_merger_dum=cr_dum(company);
own_dum=pre_merger_dum(1:7,:);



load c:/data/restat/analysis_data/postoiloutside;
meanprice=data(:,1);
meanshare=data(:,2);
meany=log(reshape(meanshare,7,10))-log(ones(7,10)-kron(ones(7,1),sum(reshape(meanshare,7,10))));


%IV has columns 1 to 9 containing the average price across each region excluding that corresponding to the current observation in each of the 9 time periods.  
%Columns 10 to 15 contain brand dummies
%Column 16 to 24 contain time dummies

%IV=[iv(:,2:2) x1(:,2:14)];
IV=[iv x1(:,4:34)];
clear iv data;
horz=['OLS    IV'];
vert=['constant  ';
      'price     ';
      'private   ';];


nmkt =27*10;     % number of markets = (# of cities)*(# of months)  %
nbrn =7;     % number of brands per market. %
cdid=kron([1:nmkt]',ones(nbrn,1));
cdindex = [nbrn:nbrn:nbrn*nmkt]';       

invA = inv([IV'*IV]);



% compute the outside good market share by market
temp = cumsum(s_jt);
sum1 = temp(cdindex,:);
sum1(2:size(sum1,1),:) = diff(sum1);
outshr = 1.0 - sum1(cdid,:);
outshr=outshr.*(outshr>=0);
y = log(s_jt) - log(outshr);

%grand means for elasticity calculations
mean_share=mean(reshape(s_jt,7,27*10)')';
mean_price=mean(reshape(x1(:,1),7,27*10)')';

%ols logit
beta_ols = inv(x1'*x1)*x1'*y;
res = y - x1*beta_ols;
%vcov_ols = 1/size(x1,1)*res'*res*inv(x1'*x1);
vcov_ols = inv(x1'*x1)*(x1'.*(ones(size(x1,2),1)*res')*(res*ones(size(x1,2),1)'.*x1))*inv(x1'*x1);
se_vcov_ols=sqrt(diag(vcov_ols));

%Newey S.E.'s
%Need to get lagged value of xe

x1_lag1=x1(1:1890-70,:);
x1_lag2=x1(1:1890-70*2,:);
x1_lag3=x1(1:1890-70*3,:);


x1_temp1=x1(70+1:1890,:);
x1_temp2=x1(70*2+1:1890,:);
x1_temp3=x1(70*3+1:1890,:);


res_lag1=res(1:1890-70,:);
res_lag2=res(1:1890-70*2,:);
res_lag3=res(1:1890-70*3,:);

res_temp1=res(70+1:1890,:);
res_temp2=res(70*2+1:1890,:);
res_temp3=res(70*3+1:1890,:);


G_0=x1'.*(ones(size(x1,2),1)*res')*(res*ones(size(x1,2),1)'.*x1);
G_1=x1_lag1'.*(ones(size(x1_lag1,2),1)*res_lag1')*(res_temp1*ones(size(x1_temp1,2),1)'.*x1_temp1);
%Need to get twice lagged value of xe
G_2=x1_lag2'.*(ones(size(x1_lag2,2),1)*res_lag2')*(res_temp2*ones(size(x1_temp2,2),1)'.*x1_temp2);
%Need to get three times lagged value of xe
G_3=x1_lag3'.*(ones(size(x1_lag3,2),1)*res_lag3')*(res_temp3*ones(size(x1_temp3,2),1)'.*x1_temp3);



vcov_ols_neweywest = inv(x1'*x1)*(G_0+.5*(G_1+G_1'))*inv(x1'*x1);
se_vcov_ols_newey=sqrt(diag(vcov_ols_neweywest));
clear x1_temp1 x1_temp2 x1_temp3 res_lag1 res_lag2 res_lag3 res_temp1 res_temp2 res_temp3 G_0 G_1 G_2 G_3;



%ols elasticities at grand means
alpha=beta_ols(1);
o1 = alpha*mean_share*mean_share';
o2 = diag(alpha*mean_share);
omega=(o2-o1);
olselas=omega.*(1./mean_share*mean_price');




% IV regressions
beta_IV = inv(x1'*IV*invA*IV'*x1)*x1'*IV*invA*IV'*y;
res = y - x1*beta_IV;
IVres = IV.*(res*ones(1,size(IV,2)));
b = IVres'*IVres;
vcov_IV = inv(x1'*IV*invA*IV'*x1)*x1'*IV*invA*b*invA*IV'*x1*inv(x1'*IV*invA*IV'*x1);
se_vcov_IV=sqrt(diag(vcov_IV));


%Newey West S.E.'s here

IV_lag1=IV(1:1890-70,:);
IV_lag2=IV(1:1890-70*2,:);
IV_lag3=IV(1:1890-70*3,:);


IV_temp1=IV(70+1:1890,:);
IV_temp2=IV(70+1:1890,:);
IV_temp3=IV(70+1:1890,:);


res_lag1=res(1:1890-70,:);
res_lag2=res(1:1890-70*2,:);
res_lag3=res(1:1890-70*3,:);

res_temp1=res(70+1:1890,:);
res_temp2=res(70+1:1890,:);
res_temp3=res(70+1:1890,:);

IVres_lag1 = IV.*(res*ones(1,size(IV,2)));
IVres_lag2 = IV.*(res*ones(1,size(IV,2)));
IVres_lag3=IV_lag3.*(res_lag3*ones(1,size(IV_lag3,2)));

G_0=IVres'*IVres;
G_1=IV_lag1'.*(ones(size(IV_lag1,2),1)*res_lag1')*(res_temp1*ones(size(IV_temp1,2),1)'.*IV_temp1);
vcov_IV_neweywest = inv(x1'*IV*invA*IV'*x1)*x1'*IV*invA*(G_0+.5*(G_1+G_1'))*invA*IV'*x1*inv(x1'*IV*invA*IV'*x1);
se_vcov_IV_neweywest=sqrt(diag(vcov_IV_neweywest));

clear IV_lag1 IV_lag2 IV_lag3 IV_temp1 IV_temp2 IV_temp3 res_lag1 res_lag2 res_lag3 res_temp1 res_temp2 res_temp3 IVres_lag1 IVres_lag2 IVres_lag3;



%iv elasticities at grand means
alpha=beta_IV(1);
o1 = alpha*mean_share*mean_share';
o2 = diag(alpha*mean_share);
omega=(o2-o1);
ivelas=omega.*(1./mean_share*mean_price');



%Approximate first-stage F-stat, leave out vcov d.f. adjustment as N is large relative to number of regresors...
% b_1st=inv(IV'*IV)*IV'*x1(:,1);
% res=x1(:,1)-IV*b_1st;
% vcov_1st=inv(IV'*IV)*(IV'.*(ones(size(IV,2),1)*res')*(res*ones(size(IV,2),1)'.*IV))*inv(IV'*IV);
% se_vcov_1st=sqrt(diag(vcov_1st));
% R=[eye(9) zeros(9,13)];
% q=9;
% F=(R*b_1st)'*inv(R*vcov_1st*R')*(R*b_1st)/q;

alpha=beta_ols(1);
pre_merger_dum=cr_dum(company);
for t=1:10
        shares=meanshare((t-1)*7+1:t*7);
    	mp=meanprice((t-1)*7+1:t*7);
        own_dum=pre_merger_dum((t-1)*7+1:t*7,:);
        o1=alpha*shares*shares';
        o2=diag(alpha*shares);
        omega((t-1)*7+1:t*7,:)=o2-o1;
        elas((t-1)*7+1:t*7,:)=omega((t-1)*7+1:t*7,:).*(1./shares*mp');
        markup((t-1)*7+1:t*7,1)=-inv(omega((t-1)*7+1:t*7,:).*(own_dum*own_dum'))*shares;
end

mc_hat=meanprice-markup;
company_post=company.*(1-(company==6))+(company==6)*4;
post_dum=cr_dum(company_post);
%post_dum=cr_dum(company);
options=optimset('Tolfun',1e-10,'TolX',1e-10,'Display','on');
first=1;
for i=1:10
	last=i*7;
    expall=exp(meany(:,i)-alpha*meanprice(first:last));
	p_star(first:last,1)=fsolve(@logitbert_eq,meanprice(first:last),options,expall,alpha,... 
    mc_hat(first:last), post_dum(first:last,:));
	first=last+1;
end
change=median(reshape(100*((p_star-meanprice)./meanprice),7,10)')




%number of bootstrap draws
sims=1000;
pchange=zeros(sims,7);
melas=zeros(7,7,sims);
for i=1:sims
    bs_beta_IV=beta_IV+chol(vcov_IV_neweywest)'*randn(size(beta_IV,1),1);
    bs_beta_ols=beta_ols+chol(vcov_ols_neweywest)'*randn(size(beta_ols,1),1);
    alpha=bs_beta_IV(1)
    markup_multi = zeros(24*10,1);
    elas = zeros(24*10,7);
    omega = zeros(24*10,7);

    % compute demand elasticities


    o1=alpha*mean_share*mean_share';
    o2=diag(alpha*mean_share);
    omega=(o2-o1);
    melas(:,:,i)=omega.*(1./mean_share*mean_price');
for t=1:10
        shares=meanshare((t-1)*7+1:t*7);
     	mp=meanprice((t-1)*7+1:t*7);
        own_dum=pre_merger_dum((t-1)*7+1:t*7,:);
        o1=alpha*shares*shares';
        o2=diag(alpha*shares);
        omega((t-1)*7+1:t*7,:)=o2-o1;
        markup((t-1)*7+1:t*7,1)=-inv(omega((t-1)*7+1:t*7,:).*(own_dum*own_dum'))*shares;
end

mc_hat=meanprice-markup;
company_post=company.*(1-(company==6))+(company==6)*4;
post_dum=cr_dum(company_post);
 options=optimset('Tolfun',1e-10,'TolX',1e-10,'Display','on');
first=1;
for j=1:10
	last=j*7;
    expall=exp(meany(:,j)-alpha*meanprice(first:last));
	p_star(first:last,1)=fsolve(@logitbert_eq,meanprice(first:last),options,expall,alpha,... 
mc_hat(first:last), post_dum(first:last,:));
	first=last+1;
end
pchange(i,:)=median(reshape(100*((p_star-meanprice)./meanprice),7,10)')
   

end
